Entropy and Correlation Functions of a Driven Quantum Spin Chain 
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We present an exact solution for a quantum spin chain driven through its critical points. Our 
approach is based on a many-body generalization of the Landau-Zener transition theory, applied 
to fermionized spin Hamiltonian. The resulting nonequilibrium state of the system, while being a 
pure quantum state, has local properties of a mixed state characterized by finite entropy density 
associated with Kibble-Zurek defects. The entropy, as well as the finite spin correlation length, are 
functions of the rate of sweep through the critical point. We analyze the anisotropic XY spin 1/2 
model evolved with a full many-body evolution operator. With the help of Toeplitz determinants 
calculus, we obtain an exact form of correlation functions. The properties of the evolved system 
undergo an abrupt change at a certain critical sweep rate, signaling formation of ordered domains. 
We link this phenomenon to the behavior of complex singularities of the Toeplitz generating function. 



I. INTRODUCTION 

Recent advances in the studies of ultracold atoms 
trapped in optical lattices have opened a new arena 
of investigation of nonequilibrium strongly correlated 
quantum systems 0, 0. These new opportunities are 
epitomized by the pioneering experiments on tunable 
Mott insulator-to-superfluid quantum phase transition, 
observed by manipulation of the optical lattice potential 
in 3d [l| and Id |3j] systems. The highly controllable 
environment and long coherence times of these systems 
provide new framework for investigation of nonequilib- 
rium dynamics of quantum critical phenomena LJ li| 

One interesting question arising in this framework has 
to do with the properties of defects produced by sweeping 
through a critical point. For the phase transitions occur- 
ring at finite temperature the defect production is de- 
scribed by Kibble-Zurek (KZ) theory 0,I3- This theory, 
which initially was applied to topological defects left be- 
hind cosmological phase transitions, and only later found 
its way in condensed matter physics, estimates the cor- 
relation length in the ordered state using a causality ar- 
gument. The correlation length serves as a measure of 
the size of the ordered domains and of typical separa- 
tion between defects. Defect production was probed in 
recent experiments employing superfluid 3 He jy, an d 
superconducting Josephson junctions pdj . 

Phase transitions in cold atom systems are character- 
ized by a high degree of coherence, which makes the dy- 
namics near the critical point essentially non-dissipative. 
The theory of defect production in this situation has to 
be modified to account for coherent dynamics. Defect 
production in quantum dynamics can be studied using 
integrable Id spin models. The Id spin models with 
varying coupling constants provide a template for many 
quantum phenomena. Realizations of such models have 
been proposed recently in Id qubit chains|l2fl and optical 
lattices 13]. The models of quantum spin quench dynam- 
ics resulting from an abrupt change of coupling constant 
which takes the system across the phase boundary, were 



considered in Refs.JfJQ]. The quench dynamics, while 
providing useful insight, do not describe the situation of 
a continuous sweep across the transition, which is ad- 
dressed in the present work. 

Besides defect production rate and density, there is an 
interesting question of the entropy associated with the 
defects. Naively, it may seem that the entropy cannot 
be produced at zero temperature by a system evolving 
unitarily in a pure state. However, if the evolved state 
is sufficiently complex, it may look entropic from a local 
point of view, i.e. if observed in a volume much smaller 
than the total system size. As we shall see, this is pre- 
cisely the case in this problem. 

In the present article we study time evolution of a 
many-body system which is swept at a constant speed 
through its quantum critical point. With the help of 
an exactly solvable Id quantum spin model with a time- 
dependent Hamiltonian we explore how the time evolu- 
tion across the critical point manifests itself in the many- 
body effects and spin correlation functions. In particu- 
lar, we analyze the relation between the sweep speed and 
spatial spin correlations, providing an extension of KZ 
scenario to the quantum critical point regime. Our an- 
alytical results are in agreement with recent numerical 
study of this problem, reported in Ref. [l^ . 

Our approach is based on a many-body generalization 
of the Landau-Zener (LZ) transition theory. In this work 
we focus on the anisotropic XY spin 1/2 chain with time- 
dependent couplings. We consider unitary evolution of 
the system, initially in the ground state, which crosses 
its equilibrium critical points. Since the Hamiltonian of 
the fermionized spin chain is quadratic, the evolution of 
the many-body state can be expressed with the help of a 
Bogoliubov transformation through a suitable set of the 
2x2 evolution problems of LZ form, one for each fermion 
momentum value. 

Our analysis reveals that the evolved system state has 
a number of interesting characteristics. Firstly, despite 
being in a pure quantum state in a global sense, its local 
properties are identical to those of a system in a mixed 
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state, characterized by finite effective temperature and 
entropy density. Although the finite entropy property of 
a pure state may seem counterintuitive, it naturally arises 
in the description of local properties, such as correlation 
functions. We shall see that the origin of finite entropy 
can be traced to coarse-graining in momentum space. On 
a more intuitive level, the system pure state can described 
as a superposition of different configurations of ordered 
domains with uniform magnetization. However, the co- 
herence of amplitudes associated with different domain 
arrangements cannot be detected locally without having 
access to the entire set of variables in the system, which 
leads to an apparent mixed state and finite entropy. 

Secondly, the transition from the adiabatic to non- 
adiabatic regime in the LZ problem, taken as a function 
of the sweep rate, depends on the momentum value of the 
fcrmionic mode. The characteristic crossover momentum 
can be associated with the inverse correlation length I 
in the KZ picture, corresponding to typical domain size. 
This approach yields a scaling relation between the cor- 
relation length and the sweep speed, I oc u -1 / 2 . This re- 
lation, obtained directly from the analysis of the many- 
body evolution operator, agrees with the KZ causality 
argument prediction. 

Lastly, due to a simple product structure of the 
evolved state, the correlation functions can be obtained 
in a closed, exact form with the help of the theory of 
Toeplitz determinants. The correlation functions exhibit 
a crossover from monotonically decreasing behavior at 
fast sweep speed, e~ r ^, to an oscillatory behavior at a 
slow speed, e~ r ^ cos(wr — ip). The oscillatory behavior, 
which appears abruptly below certain sweep speed value, 
corresponds to alternate magnetization signs in neigh- 
boring ordered domains (see Fig. The spatial period 
2?:/ to gives characteristic domain size. The parameters £, 
lo and ip exhibit a singularity at the critical sweep speed, 
which is analyzed and explained in the Toeplitz determi- 
nant framework via evolution of zeroes of the generating 
function in a complex plane. 

The plan of this article is as follows. We start with 
analyzing the full many-body evolution operator of the 
XY spin chain with the help of Jordan- Wigner fermion- 
ization and reduction to the LZ transition problem in 
each fcrmion momentum subspace (Sec. [HJ. Next, in 
Sec. IIVI we show that in a macroscopic system (number 
of sites N — * oo), a non-equilibrium steady state (NESS) 
emerges at late times. This is a mixed state characterized 
by a density matrix with finite entropy which depends 
on the sweep speed. The state of a mixed character ap- 
pears due to decoherence intrinsic to the many-body LZ 
process, without any external decoherence effects. Tech- 
nically, the mixed state arises as a result of taking the 
large N limit in the correlation functions for spins sepa- 
rated by distances much less than the system size, r -C N. 
This procedure allows to eliminate the rapidly oscillating 
terms in the correlation functions, which would disappear 
in a real system as a result of physical decoherence pro- 
cesses, even if the latter are extremely weak. The entropy 




,1 1 1 1 1 1 1 1 1 1 1 

2 4 6 8 10 12 14 16 18 20 
Distance 



FIG. 1: Spin correlation function schematic position depen- 
dence for slow sweep speed and corresponding typical arrange- 
ment of Kibble-Zurek domains. The correlation length and 
oscillation period are controlled by domain size. 



of NESS is analyzed in Sec. E\ 

The density matrix description of NESS is subse- 
quently used in Sees. IVII and IVIIII to characterize or- 
dering and analyze correlation functions. The method 
employed in analytic calculation uses some results from 
the theory of Toeplitz determinants which are reviewed 
in Appendix A. We obtain the asymptotics of equal- 
time spin correlators in the NESS which have non-trivial 
crossover behavior as a function of the sweep rate. Both 
numerical and analytical results are presented, compared, 
and found to be in agreement. 



II. SPIN CHAIN DYNAMICS 

In this section, we consider a quantum XY spin 1/2 
chain in time-dependent transverse field, described by 
the Hamiltonian 

1 N 

x=l 

where N is the number of sites. The anisotropic coupling 
values are 

Ji = J(l+ 7 )/2, J 2 =J(l-7)/2, h(t)=vt (2.1) 

Here J = A(Ji + J 2 ) is the average coupling and 7 = 
{J\ — J2)/ { J1+J2) is the anisotropy parameter. Note that 
the values 7 = 0, ±1 describe the isotropic XY model and 
the Ising model, respectively. (Without loss of generality, 
we assume J > 0.) 

In this article, the problem (|2.1|) is considered with pe- 
riodic boundary conditions, i.e. x = N + 1 is identified 
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FIG. 2: Zero-temperature phase diagram of anisotropic XY 
model adapted from Ref.jTfj. The lines of critical points, 
h — ±J and 7 = with \h/J\ < 1 are marked by dashed lines, 
the circular domain (h/J) 2 +j 2 < 1 is marked by dotted line 
(see text). The evolution trajectory of the system 12. II due 
to time-dependent h(t) is shown by solid line. 



with x = 1. Other choices, such as open boundary condi- 
tions, are possible. While the properties of interest in the 
large N limit will be insensitive to the form of boundary 
conditions, periodic boundary conditions will make the 
intermediate steps of calculations more transparent. 

The time-dependent transverse field h(t) defines the 
evolution in the equilibrium system phase space which 
starts from and ends at the state in which the external 
field h(t) is much larger than the couplings Jx,i (FigHJ. 
Thus in the asymptotic ground states at t — > ±oo the 



spins are fully polarized: ij) 

',(0) 



(0) 



(... [Hi ...) and 

tjj+'/x, = (... TTTT •••)• A fully adiabatic time evolution 
(with negligible speed v — dh/dt) would transform the 

initial state ip-lo the state V'+oo- This would also de- 
scribe physical evolution at a finite but sufficiently slow 
speed, provided that the ground and excited states are 
separated by a finite gap at all times. However, if the 
evolution takes the system through a critical point, where 
the gap vanishes, the nonadiabatic effects inevitably give 
rise to a state much more complex than (... TTTT •••)• 
To analyze the time-dependent state we evaluate the 

evolution operator Ut = Texp (-i J T T H(t)dt) , using 
Schrodinger representation. We choose a long evolution 
time interval, —T<t<T,so that 



T»/ c 



J/v, 



(2.2) 



where 2ig is the transit time between the critical lines 
h = ±J (Fig. |2J). Since the effect of the couplings J1.2 is 
important only during a relatively short time interval of 
order tQ, when h(t) ~ J1.2, one expects the results to be 
fairly insensitive to the specific value of T. Indeed, as we 
shall discover shortly, in the limit described by Ea. (|2.2|) 
universal results will arise. 

The model ()2.1f) has a long history dating back to 
the original solution of the equilibrium model by Lieb, 
Schulz, and Mattis ^} w ho obtained an exact solu- 



tion using Jordan- Wigner fermionization. Let us recall 
the basic features of the phase diagram in equilibrium. 
Barouch and McCoy 0] obtained the phase diagram by 
considering spin correlators in the ground state. These 
results were subsequently extended by Tracy and Vaidya 
[lH [Tg| and further generalized in Refs.pol I2H ] which 
employ quantum inverse scattering technique. 

For reader's convenience, here we summarize the zero- 
temperature equilibrium phase diagram as a function 
of h/J and 7 in Fig. [21 The system exhibits spontaneous 
ferromagnetic Ising order for —J<h<J, (antiferro- 
magnetic for J < 0) and can be described for \h\ > J as 
disordered, or paramagnetic. The lines of critical points 
h = ±J, separating these regimes, are in the Ising uni- 
versality class. The gap in the excitation spectrum 



e(k) =±((h+ Jcosfc) 2 +7V 2 sin 2 k) 



,1/2 



(2.3) 



vanishes on the critical lines. Outside the circular do- 
main marked in Fig. |3 j 2 + h 2 / J 2 > 1, the correlators 
in the ground state exhibit Ising-like pure exponential 
decay. In contrast, for 7 2 + h 2 /J 2 < 1 the correlators 
have oscillatory subleading terms. The ground state on 
the circle 7 2 + h 2 / J 2 = 1 is a direct product of single-site 
spin states 22] . On the 7 = line (Ji = J2) the Hamilto- 
nian is isotropic. In this case, in the interval — J < h < J 
the ground state is quantum critical. 

For our choice of the time-dependent field, the system 
is deep in the disordered phase at both the early and late 
times, \h(t ~ ±T)\ ^> J- At such times the instantaneous 
eigenstates of H (t) evolve quasi- adiabatically, with a pure 
phase factor. However, at intermediate times t ~ tQ 
we expect non-trivial dynamics as the system enters the 
phase with spontaneous Ising order, — J < h < J, passing 
through the critical points at h(t) = ± J. 

Our exact solution of the dynamical problem is a direct 
generalization of the equilibrium solution. We employ the 
time-independent Jordan- Wigner string variables 



n 



(2.4) 



x' <x 



In the Ising limit 7 = 1, the quantities t x are dual to the 
g\ and represent so-called disorder variables [2^|. With 
the help of t x we define spinless fermionic operators 

a x Tx^x ' a x = T x@x 5 

with <y~ = ^(al±ia 2 ) the raising and lowering operators. 
The fermionized Hamiltonian is quadratic: 

N 

TL = *^2A x a+a x+ i + B x a x a x +i+h-.c.-2h(t)a+ax, (2.5) 

x = l 

where we subtracted a constant Eq = Nh(t). Here the 
couplings A x — J 1 + J 2 = J , B x = J 2 — J\ — —7 J are 
the same for all 1 < x < N, and 



A, 



Jtn, B x 



-N 



-^Jtn. 



(2.6) 
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The string operator tn can be expressed as exp(wnA/), 
where Af = Ylx=i a t a x IS ^ ne total fermion number. The 
complication due to the presence of the operator-valued 
couplings (|2.6|l in the Hamiltonian 1|2.5|) turns out to be 
inessential 16]. In fact, since different terms of Ea. (|2.5(l 
either conserve the fermion number Af, or change it by 
±2, the operator tm is a constant of motion, [tjv, H] = 0. 
This allows to replace tjv by the c-number equal to its 
value in the initial state: tn = (— 1)^. Thus we obtain a 
truly quadratic translationally invariant Hamiltonian in 
the fermion representation with periodic or antipcriodic 
boundary conditions, depending on the parity of N. 

It will be convenient to write fermionic operators using 
two-component vectors, 



C k = 



a k 



ik.r t 



(2.7) 



with fc = 27rm/iV, where m is integer or half- integer, 
depending on the parity of N. The fermionizcd Hamilto- 
nian, in the momentum representation (|2.7() . splits into 
a sum of independent terms, H(t) — (^2 k>0 Hk(t)) + 
E' , where each term operates in the four-dimensional 
Hilbert space associated with the momentum states fc, 
— fc filled with different numbers of fermions, and E' = 
12k>o ^cos/c is a constant. The operators H k (t) are bi- 
linear in Cfc and have the form 



H k {t) 



fh(t) + J cos k ij J sin k 
k \ -i-fjsink -h(t) - JcoskJ V 

(2.8) 

which conserves k due to translational invariance. Also, 
H k conserves the fermion occupancy number n k = 
a k ak a ^fe a -fc U P to ±2 (i.e. the parity of n k ) sepa- 
rately within each fc-subspace (fc, — fc). 



III. MANY-BODY LANDAU-ZENER 
TRANSITION 

Using the representation (|2.8|l we can write the full 
many-body evolution operator as a tensor product of par- 
tial evolution operators acting in the (fc,— fc) subspaces: 



U(t) = 60 U k (t) , 17* (t) = Texp -i / H k (t')dt 



k>0 



—T 



(3.1) 

To obtain we consider the basis in the fc, — k subspace 
generated by the a k vacuum a k \0} = as follows: 

|0), |fc,-fe)=4a! fc |0) 
|fc)=4|0). |-fc)=aL fe |0) 

The latter two states | ± fc) of occupancy one are eigen- 
states of the Hamiltonian 12.8fl : 

H k (t)\±k) = (h(t) + Jcosfc)| ±fc). 



(This follows from conservation of fc and the parity of 
n k .) Thus each of the states | ± fc) evolves in time with 
a phase factor, | ± k){t) = e~ lv{t) \ ± fc), with 



~dt 



h(t) + J cos fc. 



(3.2) 



The other two states, |0) and |fc,— fc), evolve as super- 
position ^ k (t) = u k (t)\0) + v k (t)\k, — fc). We denote the 
corresponding 2x2 evolution operator as S k (t). 

This discussion can be summarized by writing the 4x4 
evolution operator U k in a block-diagonal form: 



S k (t) 






-i<p(t)i 



(3.3) 



with 1 a 2 x 2 identity operator. The first and the second 
block correspond to the states |0), |fc, — fc) and | ± fc), 
respectively. 

To describe S k (t), we project the Hamiltonian H k (t) 
on the subspace |0), |fc,— fc), which gives an evolution 
equation for u k (t), v k (t) as follows: 

.„ T (h(t) + Jcosfc — 2i7Jsinfc 
\ 2«7Jsmfc 



fc(t)-jco8*y** (3 - 4) 



The form of Ea. l|3.4H is identical to that of the LZ transi- 
tion problem [24L |25j for two levels evolving linearly with 
time through an avoided crossing of size = 27 J| sin fc|. 

The result of the evolution defined by Ea. l|3.4|l can 
be represented as a 2 x 2 unitary matrix which de- 
pends on the Landau-Zener adiabaticity parameter a k = 
|Afe| 2 /wi2, where V12 = 2dh/dt = 2v is the relative veloc- 
ity of the levels. The parameter a k is small for fast level 
crossing and large for slow crossing. In our case, we have 

a k = (4:-f 2 J 2 /2v) sin 2 fc = zsin 2 fc, 

where we introduced the dimensionless parameter 

z = 2 7 2 J 2 /v (3.5) 

to be used throughout the rest of the paper. 

The evolution matrix for the LZ problem can be ob- 
tained exactly in analytic form. In the limit of the total 
evolution time long compared to the level crossing time 
(realized in our case, since T/fq S> 1), one can write 
the evolution operator S k explicitly in terms of the LZ 
transiton amplitudes 



r k 



sk = -sgn(fcWl - 



(3.6) 



The long-time asymptotic form of the matrix S k (e.g., 
see Ref . [261] ) is as follows: 



S k 



r k e 
s k e 



-s k e 



i>ii,. 



rke 



where the time-dependent phases are 

<pk = fk(xt) - fk(%k) 
Vk = 
fk{x) = 



(3.7) 



/feOfc ) + fk{x k ) +tt/4 - &YgT{ia k ) 
x 2 /4 + a k \n\x\ + 0(x- 2 ) 
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FIG. 3: LZ probability 1)3.9^ of remaining in the initial state 
for z/ = 0.1, 1, 10 (from top to bottom). The dashed line 
marks pk = 0.5. Note the regions near k = 0, ±7r (critical 
modes) where LZ transition does not take place even at a 
slow sweep speed z/z* 3> 1. 



Here T(x) is the gamma function and 

x k (t) = 2(vt + Jcosk)/v 1/2 (3.8) 

is dimensionless time. Note that in the long time limit 
only the phases (fik, r\h depend on time, quickly grow- 
ing as a function of T, while the amplitudes s& be- 
come time- independent, approaching the asymptotic val- 
ues 

Since the states 0), | ± k), \k, — k) are invariant (up to 
a phase factor) at t — > ±oo, with LZ transitions between 
0) and \k, — k) happening only at times t ~ J/v, the 
asymptotic matrix Sk can be used to describe transitions 
resulting from the time evolution. In Fig. [31 we plot the 
probability 

Pk = M 2 = e- 2na " = e -2T*sin 2 fe (3 g) 

for the system, evolving from the state |0) at t = —T, to 
remain in this state at late time t = T. (The quantity 
(|3.9|) also describes the probability of the state \k, —k) to 
remain itself.) The top curve in Fig. [3] corresponds to 
small z (fast sweep rate v) when the levels cross quickly 
and the transition probability is small. The transition 
probability increases at larger z, with fully adiabatic 
regime reached for typical values of k at very large z. 
In this limit, the systems performs a nearly complete 
transfer of population from the initial state |0) to the 
state \k, — k), which in the spin language corresponds to 
spin orientation reversal — > — cr^. This behavior is 
illustrated by the lower curve in Fig. |3J In this case, 
while the majority of the modes evolve adiabatically to 
the final state \k,—k), a small fraction of the modes with 
k close to 0, ±7r evolve nonadiabatically. These modes 
remain stuck in the the initial state |0), for pf. w 1, or 
form a superposition of the states |0) and \k,—k) with 
comparable weights, for pta 1/2 (see Fig. [3J. 



To characterize the degree of adiabaticity of different 
modes, it is convenient to define a special value of z which 
will be of importance in the discussion below: 

z* = ^- = 0.110... (3.10) 

As Fig. illustrates, at z = z* the curve pk is tangent 
to the p = 1/2 line at k = ±7r/2. As we shall see in 
Scc lIVI the modes with pk = 1/2 are the ones for which 
the decoherence due to partition at LZ transition is the 
strongest. These modes at large t evolve as an equal 
weight superposition u(t)|0) + v(t)\k, — k) with \u(t)\ — 
\v(t)\ and relative phase rapidly changing in time. The 
oscillatory phase factors will be identified below with the 
source of intrinsic decoherence. 

In addition, we shall see in Secs lVIl and IVIIll that the 
value z = z*, which marks the appearance of the modes 
with pk = 1/2, is also special in another way. We shall 
find that the spin correlation functions in the final state 
undergo an abrupt change at the sweep speed value cor- 
responding to z = z*, from monotonic at z < z* to os- 
cillatory at z > z*. Interestingly, this transition in the 
correlation function behavior occurs at the same speed 
value which corresponds to the largest phase space of the 
modes with pk = 1/2. 

At fixed z, the degree of adiabaticity for a particular 
mode is quite sensitive to the value of k. As illustrated 
in Fig. |31 due to the sin 2 k dependence in pk, the adia- 
batic regime for the modes with different k is reached at 
different values of the sweep speed, z sin 2 k 3> z*. In par- 
ticular, for the modes with k sufficiently close to and 
±7r the transition is adiabatic only at very large z. These 
modes are special since they are gapless on the critical 
lines h = ± J of the equilibrium phase diagram, crossed 
by the evolution trajectory (Fig. |2J). Such critical modes, 
characterized by small excitation frequency, vanishing at 
k = 0, ±7r, are not able to react to field sweep with fi- 
nite velocity v, no matter how small the latter is. For the 
whole system, the nonadiabatic behavior of the k = 0, ±7r 
modes means that the spin reversal is incomplete even at 
very slow sweep. The fraction of the spins that do not 
accomplish reversal, at large z can be estimated as 

An = J2 P k ~ \ I e ~ 2 ™ k * dk = (^ 2 z)- 1/2 . (3.11) 

The density of defects An has an inverse square root de- 
pendence on the sweep speed v. By order of magnitude, 
the estimate (I3.11f) can be obtained also from the momen- 
tum value k ~ (z^/z) 1 / 2 corresponding to the crossover 
at p k ~ 1/2. 

Our result H3.11(l for An can be compared to the es- 
timate following from the KZ causality argument 0, 0] , 
which predicts the domains of the ordered phase of size 

£ = ct (3.12) 

where c is the velocity of gapless excitations at the critical 
point and r is the characteristic transit time. In our case, 
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from the excitation spectrum (|2.3|l . at the critical points 
h = ± J the velocity is c = 7J. The transit time for the 
/c-mode can be estimated as the time of sweeping across 
the gap: r k ~ A^/v, where A& = ck. After identifying 
I with 1/k, Ea. pil|l becomes i = c 2 /(v£), yielding I vs. 
v dependence 

e = c/y/v- (3.13) 

The —1/2 power law scaling is in agreement with the 
result (|3.11() , which confirms the KZ scenario 0, § for 
Id spin chain and links it to the many-body LZ transi- 
tion. Similar observations were made in a recent numer- 
ical study of a spin model in a finite size system |15|. 



IV. DECOHERENCE DUE TO TRANSIT 
THROUGH CRITICAL POINT 

Here we discuss the phenomenon of intrinsic decoher- 
ence resulting from massive production of spin excita- 
tions at a sweep through critical point. We start with 
noting that the evolution during — T < t < T, taken 
formally, is manifestly unitary and preserves all phase 
relationships. For the density matrix of the entire sys- 
tem, the evolution ip = [p, H] starting with a pure state 
p(t = —T) = |0iv) (Oat I of A spins obtains a pure state: 

p(N,T) = p(t = T) = U T \0 N )(0 N \tf T (4.1) 

However, we shall see that some of the phases in the 
density matrix (|4.1(l develop rapid oscillation at large T. 
The phase growing with T will be found to depend on the 
momentum k so that the oscillatory part of p averages to 
zero after integration over k in all local correlators. 

It is beneficial to identify the oscillatory terms and 
suppress them early in calculation. This can be achieved 
by replacing the unitarily evolved state l|4. 1J) by a deco- 
hered state from which the rapidly oscillating terms are 
removed. This procedure, which will be seen to describe 
correctly the local properties of the evolved system, leads 
to the notion of non-equilibrium steady state (NESS). Be- 
sides the benefit of simplicity, early appearance of NESS 
in the analysis also helps to develop intuition about how 
the results of evolution depend on various parameters, 
such as the sweep rate and coupling strength. 

Alternatively, one could proceed more formally, car- 
rying the oscillatory terms in p over and then arguing 
that they drop out in the limit of long time T and large 
systems size. For that, one would have to include the 
effect of some auxiliary phyical decoherence mechanism 
and obtain suppression of the oscillatory terms indepen- 
dent of the strength of decoherence effect, no matter how 
weak the latter is. Instead, we choose to build the NESS 
and its decohered density matrix prior to analyzing the 
correlators. 

We shall focus on the observables, i.e. spin correla- 
tors, which are more physical quantities than the full 



many-body density matrix. Let us consider correlators 
in position space within a block [1, n\: 

(A{x)...A'(x')}, 1 < x,...,x' < n, (4.2) 

where A...A' are local observables given by products of a 
finite number of fermion operators. In the discussion be- 
low the intermediate length scale n will be much smaller 
than the system size, n <C A. These correlators can be 
evaluated with the help of a reduced density matrix ob- 
tained by tracing out all spin variables outside the block 
1 < x < n. The resulting density matrix describes only 
the 2™ spin states in the block: 

p(n,T) = Tr N _ n [p{ A, T)\ (4.3) 

where Trjv_„ denotes integration of the N — n spins out- 
side the block 1 < x < n, and T is evolution time. The 
reduced density matrix adequately describes the correla- 
tors and other properties at distances shorter than n. 

Next, we consider how by taking the three scales A, 
T, n to infinity in proper order one arrives at the NESS. 
First, we take the thermodynamic limit A — * 00 to elim- 
inate recurrence times of order of level spacing for finite- 
size systems. Second, we take the long time limit T — > 00 
to suppress oscillations and arrive at a steady state. Fi- 
nally, we take the long-wavelength limit n — > 00 and 
obtain the decohered density matrix 

po — lim lim lim prj(N,T,n) (4-4) 

n^oo T—>oo N-^oc 

which describes NESS in an infinite system. 

Not all phase relationships of the pure state density 
matrix, Eq. (|4.1|) . survive the limiting procedure (|4.4|l . 
We describe this process as decoherence by analogy with 
loss of phase information for density matrices of open 
quantum systems coupled explicitly to an environment 
[271 128| . In contrast to the latter, however, the decoher- 
ence described by l|4.4|l is of an intrinsic origin, arising 
from the spin chain acting as both the system undergo- 
ing decoherence and the environment that induces it. An 
implicit separation between the two emerges only when 
considering correlators in contrast to the explicit separa- 
tion in open quantum systems. 

We shall use the fermion representation constructed 
above to evaluate po- Formally, this restricts our the- 
ory to the correlators of the form (|4.2I) with the observ- 
ables A. ..A' all taken at equal times. In the fermion 
representation, the full density matrix p(N, T) decou- 
ples into a tensor product over the k,— k subspaces: 
p(N,T) = <S> k>Q Rk with a 4 x 4 matrix R k . The lat- 
ter has nonzero elements only between the states |0) and 
\k, — k), since the amplitude of the states | ± fc), which 
is zero in the initial state, cannot change with time [see 
Ea. (|3.3|) ]. Thus within each k — k subspace the density 
matrix R k is effectively 2x2, restricted to the subspace 
|0), \k, — k) where it is nonzero: 

,(A,T)=(g)p fel , = ( p ;r-t) ^ 

fe>0 v 7 



7 



where p k is evaluated as 5 , fc|0)(0|S , fc 1 . Here p k is given 
by Eq.JsI33, and 

qk = r kSk e l{ ^+^ (4.6) 

are obtained from the S-matrix 1)3. 7|l . 

Now, let us consider correlation functions in the 
fcrmion representation. Since the Hamiltonian is 
quadratic in this representation, the state p(N,T), ob- 
tained by evolution of the t = —T fermion vacuum, is 
of a gaussian form. This allows to employ Wick's theo- 
rem to write any correlator as a sum of products of pair 
correlators. Thus an arbitrary local observable can be 
expressed in terms of the 2x2 matrix of pair correlators 

G(x, x',N, T) = (C x Cl) = Tr[p(7V, T)C x Cl] (4.7) 

while (C x ) = for a gaussian fermion state. 

Using Ea. (|4.7|) we can obtain the decohered matrix po 
by demanding that it reproduces G(x, x' , N, T) under the 
limits in Ea. (|4.4|) . Taking N — > oo first, we write the 
result as an integral over a continuous k variable: 

G (^', TC ,r)=£| e -*-.-)(; ;i !g (4.8, 

Turning to the T — > oo limit, we note that, while p k 
and the modulus \q k \ approach the asymptotic LZ values 
exponentially quickly, the phase of q k exhibits oscillations 
as a function of time T. To the leading order, at large T 
we have 

+ Vk ~ vT 2 + 2 JTcos k + O(lnT). (4.9) 

It is crucial that this phase has a cosfc dependence on k. 
Due to the fc-dependent phase factor e^ ipk+rik \ with the 
oscillations becoming very fast at large T, the integral of 
qk over k in Eq. 1)4.8(1 vanishes in the limit T — > oo (which 
practically means T/tQ 3> 1). This argument shows that 
the off-diagonal elements of the correlator G(x, x', oo, oo) 
vanish for arbitrary x, x'. The result 

G (^»,»)=/"je— f o » ,_g 

(4.10) 

means that can simply ignore the oscillatory terms by 
setting qk = in all correlation functions. We point 
out that such a result agrees with the intuition that the 
rapidly oscillating off-diagonal matrix elements of p van- 
ish due to arbitrarily small decoherence, and thus can be 
ignored in all correlation functions. 

Applying the q k = rule to p(N,T) given by Ea. (|4.5|) 
we obtain the decohered density matrix po as a product 
of diagonal 2x2 matrices, restricted to the subspace |0), 
\k,-k): 

PD=<S> PD,k, PD,k = ( P q 1 _ J (4.11) 

fc>0 ^ ' 



With such identification, the decohered pair correlator 
is G(x, x', oo, oo) = Ti-[pdC x C x ,], as required. The 
relation of the higher order correlators with the pair 
correlators via Wick's theorem decomposition, including 
fermionic signs, remains unchanged. 

Finally, we note that p k = e~ 27r2sm k as a function 
of k has periodicity 7r, while q k has periodicity 2tt (see 
Eas. l|3.6) ). H3.7|l '). Thus the correlation functions of the 
decohered state po, obtained by setting q k = 0, acquire 
even/odd sublattice structure in position space. The cor- 
relators 14.10(1 vanish if x and x' belong to different sub- 
lattices: 

x-x' = 2n+l: G(x,x')=0 (4.12) 
where 

p k e- ikn ^=e-™I n (irz) (4.13) 

-7T 271 

with I n (x) the modified Bessel function of the first kind. 
The decoupling of the even and odd sublattice in the 
decohered state, manifest in Ea. ((4.12(l . indicates that the 
decohered density matrix factorizes as 

Pd=Pe®Po, (4-14) 

where each pg, po acts only on the even (odd) sublattice. 
This factorization will be used below in the analysis of 
spin correlation functions. 

V. ENTROPY OF THE DECOHERED STATE 

The necessity of transition from the pure state to 
NESS, characterized by the decohered density matrix po, 
can be inferred without reference to pair correlators, by 
employing the procedure of coarse-graining in momen- 
tum space. Let us consider the evolved pure state den- 
sity matrix p(N, T), Eg. ((4. 5(1 . While the diagonal matrix 
elements of p(N, T) are smooth functions of k and in- 
dependent of T, the off-diagonal elements between |0) 
and \k,—k) rapidly oscillate as functions of both k and 
T. The oscillation fc-dependence, described by the phase 
factors e ± 2lJTcosk [ see Ea. (l4.9ll ]. becomes increasingly 
more rapid at large T. This property makes the oscilla- 
tory terms very sensitive to coarse-graining in k space: 
They vanish after intergrating over any small interval 
Afc <C 1 which is large compared to (JT) _1 . This argu- 
ment, applied above to individual correlators evaluated 
at finite separation in real space using the integral repre- 
sentation 1(4.10(1 . can also be applied to the entire density 
matrix. The coarse-graining selects the matrix elements 
of p(N, T) which are smooth in k, suppressing the oscil- 
lating parts. Only the diagonal elements of p(N,T) sur- 
vive in po, consistent with the interpretation that the su- 
perpositions of the |0) and |fc, — k) states decohere into a 
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statistical mixture. Using the language of open quantum 
systems [28|, one can identify the instantaneous eigen- 
states |0), | ft, — ft) of H(t — > +00) with the pointer states 
which survive decoherence. 

To quantify the amount of information lost in the de- 
coherence process [28J, we consider the von Neumann 
entropy of the system, S = — tr pu lnp^. [It will be 
more convenient to use natural base In instead of a more 
standard log 2 .] An expression for the entropy density 
s = S/N follows directly from the form (14. lip of po- 



(Pk Inpfc + (1 - Pk) ln(l - p k )) 



dk 
2^ 



(5.1) 



Using Taylor series for ln(l — p^) and evaluating the in- 
tegral for each term, we obtain 



s = (ttz + l)po(z) - nzpi(z) - 



p (z(m + 1)) 
m(m + 1) 



(5.2) 



with p n {z) given by Eq. l|4.13|) . The entropy (|5.2[1 as a 
function of sweep rate is plotted in Fig. We note that 
s tends to zero in the limit of small and large z, since 
for such z the dynamics gives rise to few superposition 
states. The function s{z) peaks near z ~ z*. 

Let us consider the limit of slow sweep speed, z»z,. 
In this case, Eci. (|5.2|) gives 



1.00 



0.75 



w 0.50 



0.25 



0.00 



0.01 




100 



FIG. 4: Entropy density s, Ea. l|5.2^ . as a function of z/z*, 
inverse sweep speed. Note that s peaks near z, and tends to 
zero for small and large z. 



than in an ideal gas. Another manifestation of partial 
ordering and correlations of the defects will be discussed 
in Secs lVll and lVlllI where we shall see that at slow sweep 
speeds the two-point spin correlation function exhibits 
spatial oscillations, with abrupt onset at z = z*. As 
illustrated by Fig. such oscillations result from quasi- 
regular arrangement of KZ domains. 



— ' ni(m + l) 3 / 2 



An w 0.0761 An, (5.3) 



where An = (27r 2 z) -1 ' 2 is the density of defects in the 
spin-reversed state (|3.11|1 . which describes the fraction 
of the spins remaining not reversed after slow evolution. 
For fast sweeps, z <C z*, using the expansion pj, = 1 — 
27rzsin 2 ft, we obtain 



s = — I z 

' — 7T 



sin 2 ft In 



^-)dfc«p ln— (5.4) 

27rz sin k J pa 



where po — ttz is the density of the defects evaluated as 
(1 — pk)dk/2TT, and c ~ 1 is a constant. In this case, 
Po describes the number of reversed spins, which is small 
at a fast sweep. 

It is interesting to compare these results to the entropy 
of a classical gas of a low density An -C 1, 



"lias 



= -An In An - (1 - An)ln(l - An) « An In— . 

An 



This agrees with the result for the fast sweep, Eq.JQJ, 
upon identification of An with po- In contrast, the value 
s obtained for slow sweep, Ea. 1)5.3(1 is small compared to 
Sgas for the same density: 

Sgas/s ~ 13.966 ln(e/An) > 1. 

Small entropy indictes that the arrangement of defects in 
the quantum system after a slow sweep is more orderly 



VI. SPIN CORRELATORS AND TOEPLITZ 
DETERMINANTS 

Here we consider the correlation functions of spin vari- 
ables er", and of the string variable t x , Ea. (|2.4l) . used 
in the fermionization transformation. We obtain exact 
expressions for these correlators in the form of Toeplitz 
determinants, which will allow to analyze them at large 
spatial separation. We shall see that the asymptotic 
behavior of the correlation functions is sensitive to the 
sweep speed, changing abruptly from a pure exponential 
decay at z < z, to an oscillatory dependence at z > z*. 
In this section, we focus on the non-trivial behavior for 
the correlators of transverse spin a\ , er 2 . , and of t x , and 
give a simple mathematical and physical picture. The de- 
tailed derivation and additional results on a x correlators 
can be found in Sec. I Villi 

It is convenient to write the quantities of interest as 
products of Majorana fermion operators A x = a x + a x , 
B x = a x — a x . (For convenience, we omit the factors 
1/V2, i/V2 often appearing in the definition of these 
operators.) The Majorana operators A Xl B x satisfy the 
algebra 

Al=A x , [A x ,Ay]+ = 5 xy (6.1) 
Bl - —B x , [B x ,B y ]+ = -5 xy (6.2) 
[A x ,B y ]+ = (6.3) 

In the fermion representation, the pair products of the 
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spin variables a x as well as the string variables r x , ap- 
pearing in the correlators, can be expressed as products 
of Majorana operators as follows: 



_l_i 

a x a x+n 



B x A x+ iB x+ i . . . A x+n _iB x+n _iA x+ {pA) 

fJ x <T K+?i = A x A x +iB x +i . . . A x + n -iB x + n -iB x +y3.5) 



^~ x ^~ x -\-n — A x B x A x +\B x +i . . . A x j rn B x j rn , 



(6.6) 



To obtain expectation values, we average the products 
of a finite number of the operators A x , B x using Wick's 
theorem and the decohered density matrix pn, Ea. H4.lljl . 
introduced in Sec lIVI 

An additional simplification occurs due to decoupling 
of the fcrmionic correlators, evaluated with the decohered 
density matrix po, into a product of separate contribu- 
tions of the even and odd sublattice, Ea. lj4.14jl . Let us 
explore this factorization for the correlator (cr x a x+2 n)- 
By regrouping the operators A x , B x , separating the parts 
corresponding to the two sublattices, we write 



<7 X <T X +2n 



(B x A x+ 2B x+ 2 ■ ■ ■ A x _)_2 n _2-Bx+2n-2-B x +2 n ) 
{A x +iB x +i . . . A x J r 2n-lB x J r 2n-l)- (6-7) 



Comparing the two expressions in parentheses to 
Eas. lj6.5j ). lj6.6l) . we see that the spin operator pair prod- 
uct (J x al +2n evaluated on the full lattice is a product of 
analogous operators o\.a\ +n and T x r x + n , each evaluated 
on a sublattice. This leads to factorization for the ex- 
pectation values since fermionic pair correlators do not 
mix different sublattices. The result can be symbolically 
written as 

(<rl<rl+2n) = (Wl^l+n)){{ T xT-x+n}) (6.8) 

where the brackets (...) describe expectation values of 
operators on the full lattice, while ((...)) refer to an ex- 
pectation value on a sublattice. Similar reasoning for 
other correlators leads to 



2 2 
J x a x+2n 



) = {{<<+n))((T x T x+n )) 



(6.9) 



(T x T x+2n ) = ((TxTx+n)) ((TxTx+n)) (6.10) 

where single (double) brackets refer to correlators on the 
full lattice (sublattice). This allows us to focus just on 
the sublattice correlators. 

With the help of fermionization, the sublattice corre- 
lators at separation n can be written in terms of n X n 
determinants of Toeplitz matrices, defined by a set of 
constant diagonals: 



D n [f] =det 



/o f-i 
/i /o 

/n-1 /n-2 



/-(n-1) 
/-(n-2) 

/o 



(6.11) 



The structure of the matrix, completely specified by the 
set of numbers /„, can be encoded in a generating func- 
tion 

/(£) = E /« = / yr n /(0 (6-12) 

n J C S 



with the contour C being the unit circle |£| = 1. The 
properties of Toeplitz determinants depend on the com- 
bination of poles, zeros and other singularities of /(£) in 
the complex plane p9| . 

In our case, the Toeplitz matrix representation is ob- 
tained by evaluating the sublattice correlators in (j6.8jl . 
(I6.9jl , (j6.10jl using fermion representation. With the help 
of Wick's theorem, all correlators can be expressed as 
polynomials of pair correlators of Majorana fermions. 
Due to the sublattice structure of pu, the nonzero pair 
averages are all of the form (a x a x >) with x — x' even. In 
addition, the expectation values (A x A x i) and (B X B X >), 
with x 7^ x', are zero due to Majorana fermion algebra. 
Only the pairs of operators B x A x i give nonzero expecta- 
tion values: 



(B X A X ,) = / e^ x - x '\l-2 Pk ) 



dk 
2^ 



(6.13) 



where p k is the LZ probability, Ea. (j3.9|) . [We note that, 
while A x A x i = B X B X ' = 1 at x = x', such combinations 
do not arise in the fermionic representation of spin vari- 
ables.] Summing over all pair contractions with appropri- 
ate fermionic signs brings the sublattice spin correlators 
to the Toeplitz determinant form: 



(( T x T x+n)) 



D n [g +hz ] 
D n \g- l > z ] 



(6.14) 
(6.15) 
(6.16) 



where the generating functions g m z are defined as 

ff ™.*(0 = -(-0 m (l-2p fc ), C = e 2lfe . (6.17) 

This form of the generating function, and, in particular, 
the origin of the factors £, can be understood as 

follows. The string of A X B X operators appearing in the 
correlator has an additional B x at the beginning and 
A x+n at the end compared to a similar string for the r 
correlator. This results in a shift of the matrix elements 
g n — * g n +i in the determinants for r compared to the one 
for a x , which translates to the mapping g(£) — > of 
the generating functions. Similar reasoning accounts for 
the factor £ -1 for the a x correlator generating function. 
The factor of 2 in the relation £ = e 2ik arises because the 
correlators are restricted to a sublattice, which makes the 
/c-dependence ir periodic rather than 2tt periodic. The 
factor — (— l) m ensures correct fermionic sign. 

The Toeplitz matrix representation allows to study the 
correlation functions numerically, since evaluating deter- 
minants on a computer is a low cost operation. However, 
as we show below, the problem can also be handled an- 
alytically. The benefit of the analytic treatment is that 
it provides a very clear and complete description of the 
behavior of the correlation functions at different sweep 
speeds, including the transition at z = z*. 
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VII. SPIN CORRELATORS ASYMPTOTICS 



We are primarily interested in the behavior of the sub- 
lattice correlators at large separation which maps to the 
large- n asymptotics of Toeplitz determinants. It is in- 
structive to recall the Szego limit theorem result for the 
Toeplitz determinant f6.11[) asymptotic behavior: 



(7.1) 



which holds when the generating function /(£) has a zero 
winding number and no singularities on the unit circle. 
The origin of the asymptotic (|7.1|l can be seen by noting 
that in this case the matrix elements f n rapidly decrease 
with \n\, and the Toeplitz matrix can be approximated 
by a band matrix. Then the result (|7.1|) naturally follows 
after closing the interval 1 < x < N into a circle and go- 
ing to Fourier representation. The question of how the 
asymptotic <|7.1|) is modified in the cases when the wind- 
ing numbers are nonzero and/or the generating function 
has singularities on the unit circle has been a subject of 
many publications. Not trying to review all the litera- 
ture, in the discussion below we will cite the available 
results, either conjectured or proven, as appropriate. 

We shall start with the simplest situation, for which 
Szego limit theorem provides a suitable framework. Let 
us consider the Toeplitz determinant representation for 
the correlator l|6.16[l with the generating function /(£) = 
g°' z (£). This function is real for |£| = 1, and thus has 
zero winding number. In this case, Ea. (|7.1|) yields 



D n [g°' z ]^e an , a. 



In 1 - 2e 



-2-Kzsin k 



dk 



The expression for a is analytic at z < z», has a singu- 
larity at z — z„, and becomes ill-defined at z > z*. To 
clarify the origin of this behavior, let us inspect zeros of 
g°' z . There is an infinite number of zeros £ = A p , A~ 1 of 
multiplicity one, with p an integer, which can be found 
from the representation 



z{l-2z*/z-x) _ i 



where x — (£ + £ 1 )/2. We obtain 



exp 



acosh 1 



In 2 

TTZ 



2ip 

z 



(7.2) 



(7.3) 



where we choose the branch of acosh(x) with positive real 
part so that |A P | < 1. Note that |A P | > \\' p \ for \p\ < \p'\, 
so that the zeros closest to the unit circle are Ao and Ag 1 , 
which satisfy 



-(A + A- 1 ) = l 



2z*/z. 



(7.4) 



The X(z) dependence has a square root singularity at z = 
z* . To specify the analyticity branch near the singularity, 




FIG. 5: Motion of the roots Ao, Aq" 1 as a function of z from 
the negative real axis for z < z» to the unit circle for z > z*. 
The direction of the arrows indicate increasing z. 



we take p = +0, with an infinitesimal positive part, in 
Ea.lTrSl). 

The z dependence of the roots (|7.4|) is illustrated in 
Fig0 Both roots are real and negative at z < z„: Ao < 
— 1 < Aq" 1 < 0. As z tends to z*, the roots move along 
the real axis towards £ = — 1, approach one another and 
merge at z = z* , then split and remain on the unit circle 
at z > z*, with Aq 1 = Aq. This leads to a singularity 
of the determinants D n [g°' z ], and thus of the correlation 
functions, at z — z*. 

To better understand the behavior near z — z*, it is 
instructive to try isolate the effect of the roots Ao, Aq . 
For that, we consider simplified generating functions 

/ (m) (0 - -(-CrAo 1 (1 - A £) (1 - Aor 1 ) (7.5) 

where m = 0, ±1 and Ao, A^ 1 are defined by Eg. 17.411 . 
The simplified expressions (|7.5|l capture most of the non- 
trivial behavior of the sublattice correlators arising at 
z w z* . Each of the functions /( m ) has only three nonzero 
Fourier coefficients /| m ' , and thus the Toeplitz matrix in 
this case is three-diagonal. One can easily calculate the 
corresponding Toeplitz determinants, obtaining 



D n [f {±1) ] = (-iy 



Dn[f {0) ] = (-1)"^- 



in+1 



- A 



-(n+1) 



Aq — A n 



(7.6) 
(7.7) 



These quantities, obtained from the simplified generating 
functions, Ea. f7.5fl . describe the qualitative behavior of 
the sublattice correlators for a\ , <j\ , and t x , according to 

Eqs . inni) , JOSJ , jHjg ■ 

The expressions (17.6(1 are independent of Ao, indicating 
a smooth behavior of a^, a\ sublattice correlators with 
z which will persist upon including the full generating 
function. The m = determinant, Ea . (f 7. 7fl . is analytic 
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as a function of Ao even at Ao = Aq~ . More interestingly, 
and somewhat unexpectedly, it is analytic as a function 
of z at z = z», since the right hand side of Ea. (|7.7(l is 
polynomial in Ao + Aq 1 . As a function of n, the expres- 
sion H7.7f) exhibits a crossover from exponential behavior 
at z < 2* to oscillatory behavior as z > z*. In addi- 
tion, it grows linearly with n exactly at z = z*. This 
crossover behavior, as well as nonanaliticity in z, persist 
upon including the full generating function. 

For comparison, let us consider the asymptotics for 
sublattice correlators obtained from the full generating 
function, as discussed in Sec. IVIIII 



««+„» « E 1 {-G) n , (a = 1,2) 



((TxT x+n )) «^l(-G) ? 



A™ +1 £ 2 



An — A n 



(7.8) 
(7.9) 



where G and E ly2 , given by Eas. H8.4[) . (|8.5l) . (|8.6|) . have a 
smooth z dependence. We note similarity of the behav- 
ior of these expressions to Ecis. l|7.6|) . H7.7|) at ~ Aq~" . 
We see that the origin of the crossover behavior in the 
sublattice correlators, resulting in nonanaliticity in z, is 
indeed the motion of the zeros Ao from the real axis to 
the unit circle. 

It will be useful to also write the sublattice correlators 
in the canonical form 

((«+„» ^A a e- n ^ cos™ (7.10) 
((r x T x+n )) « A r e~ n ^ T cos(Lu T n - tp T ) (7.11) 

(q = 1,2). The parameters appearing in these expres- 
sions, the amplitudes A aT , the correlation lengths £ a ,r, 
the wavenumber of spatial oscillations u> T , and the phase 
ip T , are plotted as a function of zj 'z* in Fig. 

The correlation lengths l a and £ T both become large 
at a slow sweep speed. At a fast sweep, the a^' 2 correla- 
tors become short-ranged, while the t x correlator is long- 
ranged. The oscillatory behavior of the t x correlator ap- 
pears abruptly at z = z*, with the spatial frequency and 
other parameters displayed in Fig. [^having non-analytic 
behavior. The character of this singularity is similar to 
that exhibited by the simplified model discussed above, 
Eas. l|7.6) ). l|7.7[) . which is controlled by the zeroes of the 
generating function nearest to the unit circle. 

Although the sublattice correlators are mathematically 
convenient, the physical content of our results becomes 
more transparent in the full lattice correlators. From the 
factorization relation, Ea. (|6.8() . since cos 7m = (— )", the 
a 1,2 correlators are simply given by 



over 



x u X +2n 



) w A a A T e~ n/t cos((7r - uj T )n + <p T ) (7.12) 



(a = 1, 2), where l~ x = l~ x + £-\ 

Now, let us discuss the physical regimes described by 
these correlations functions. In the time evolution con- 
sidered here, the system is driven from the disordered 
through the ordered phase and back into the disordered 
phase. In equilibrium, the correlations of a\: 2 are ab- 
sent at the early and late times, i.e. in the disordered 
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FIG. 6: The sublattice correlation function parameters, 
Eas, H7.10M7.11fl . dependence on the inverse sweep speed 
z/z*: (a) the correlation lengths l a% £ T ; (b) the frequency 
uj t and phase shift tp T ; (c) the amplitudes A a , A T . Shown 
are the analytical dependences obtained from Eqs. 17.8ll . 17.9L 
which were verified by evaluating Toeplitz determinants nu- 
merically. 



phases, but can build up at intermediate times when the 
system is in the ordered phase. The simplest situation 
arises at small z, i.e. high sweep speed. In this case, 
all the modes in the system can be treated in a sudden 
approximation. There is very little time for correlations 
in the ordered phase to build up, which results in very 
short range correlations described by exponential decay 
with a small correlation length. 

In contrast, large z describes the slow sweep speed 
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regime, when the dynamics becomes more adiabatic. 
However, full adiabaticity cannot be reached for a system 
driven across quantum critical points where the gap van- 
ishes. The build-up of correlations upon crossing the first 
quantum critical point from the disordered to ordered 
phase, h = —J, can be understood in the KZ frame- 
work as appearance of ordered domains of size I « cj\Jv, 
Eo. (|3.13fl . The length scale I characterizes separation be- 
tween defects of the ordered state, resulting from nona- 
diabaticity at crossing the critical point. The defects for 
the ferromagnetically ordered state, describing our sys- 
tem in equilibrium at — J< h(t) < J, are domain walls 
separating domains with opposite magnetization. The 
magnetization sign alternation in the domains leads to 
an oscillatory behavior of the correlators on top of ex- 
ponential decay, as illustrated in Fig. ^ Subsequent 
crossing into the disordered phase at h = J then leads 
to suppression of the correlations built up in the ordered 
phase. This is consistent with the behavior of the full lat- 
tice correlators where both the correlation length £ and 
the spatial period 2-k/(tt — lu t ) grow as z 1 / 2 at large z, 
while the correlator amplitude A r A a goes to zero. 

The crossover between the two regimes, occurring at 
z = z. t , corresponds to sweep speed of the order of the in- 
verse bandwidth. As in the discussion of Eas. lf7.6f) . (f7.7l) . 
the behavior of the correlation functions near z — z* at 
fixed n is analytic and is described as a smooth crossover. 
The apparently singular behavior in Fig. [S]is analogous 
to Stoke's phenomenon |30| for asymptotic series, where 
the coefficients of an asymptotic expansion of a function 
may not be analytic in some parameters even when the 
function itself is analytic in those parameters. 

VIII. SPIN CORRELATORS II 

In this section we outline the details of derivation of 
the results discussed above. A general procedure for cal- 
culating the asymptotics for Toeplitz determinants from 
the structure of the singularities in the generating func- 
tion is described in Appendix IA1 We note that, while this 
procedure in its most general form is only a conjecture, 
it is a reasonable extension of known rigorous results. 
Moreover, since our generating function, Eq. H6.17fl . has 



only simple zeros of integer order, the approach used be- 
low stands on a firm ground: In this discussed in 
Appendix [XJ our procedure follows from a rigorous result 
of Ref. 31]. In addition, we have compared our analytic 
results to the correlation functions obtained from direct 
numerical evaluation of the Toeplitz determinants, and 
found them to be in full agreement. 

As discussed above, among all roots of our generating 
function, Ea. (|7.3f) . one pair, Ao and Aq , plays a special 
role. We write the generating function g m ' z in the form 

9 m ' z (0 = -(-0 m A 1 (1 - AqO (1 - Aor 1 ) e h «> (8.1) 

which isolates these most relevant roots into a factor 
identical to the simplified generating function discussed 
above, Eo. H7.5fl . The remaining part, e^vO , has the form 

e H0 = iwfe+r 1 )/* nz TT z(l - XpQjl - x^- 1 ) 
V2^ P Vo 4|p|A p 

(8.2) 

It is explicit in this expression that e h ^ has all its sin- 
gularities located further away from the unit circle than 
Ao and Ag -1 . The expression for h(£) can be written in a 
more compact form: 

A «» = 1 °( 2(1-2*,/*-*) ) (8 ' 3) 

where x = (£ + £ _1 )/2. 

We obtain the correlator asymptotics given in 
Eqs. Q , (Z!13l either by using the result of Ref. or 
by the more general method of Appendix^] The latter 
procedure involves a contour C that passes through the 
two roots Ao, Aq" 1 closest to the unit circle. (This con- 
tour does not have to be a circle when |Ao| ^ 1.) We 
isolate the contributions of Ao and Ag -1 , and incorporate 
the rest of the generating function into a part smooth off 
of C, denoted by h(£). In the contribution of to the 
quantities in the asymptotics, given by contour integrals 
over C, we can deform C to the unit circle. Finally, we 
reparameterize the complex variable £ on the unit circle 
in the contour integral with x = cos 8, £ = e 10 . This 
yields expressions for the parameters G, Ei of the form 



\uG = /,,= / ^-^t=h(x) (8.4) 

l r i 



x-y V 1 - y 



x > — z—. — \/i — — ( 8 - 5 ) 



ln£ 2 = £ M AF-A -») = -ln(^)e^-,)-i/ 1 ~^=Kx) ^ ^ (8-6) 



where and €)(x) is the step function. Here h n are the Fourier coefficients of h(£) = X)n' l «£™> an d the function 
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h(x) is just h(£) reparameterized with x — (£ + £ _1 )/2. 
The integrals are in the principal value sense where ap- 
propriate. Both G and Ei are positive and smooth in 
z. However, E% is real for z < z* and a pure phase for 
z > Zit which is the same behavior as Ao- 

In passing from correlators in the form of 
Eqs.C3l,CiS) to the canonical form of Eqs. jZHS^CHH 
we have to formally drop the Aq for z < z* since it is 
subleading compared to A^ n . This gives the following 
relation between the parameters: 



e- 1 

lnA T 



Ei 

-hxG 
Re In 



■7) 



Ei 



-Re In 



(8.9) 
(8.10) 



G_ 

V 

uj t = (tt - Im In A )6(z - z«) (8.11) 
ip r = (Im \n{\ Q E 2 ) - 7r/2)0(z- z*) (8.12) 

where a = 1,2 and 0(a;) is the step function. We note 
that tv T and t/j r are nonzero only at z > z*. Also, the 
two correlation lengths t a and £ T are equal at z > z* and 
differ at z < z» (see FigEJl. 

Now, let us consider the behavior at large z, describ- 
ing slow sweep speed. As we noted earlier, in this case 
the length scales £ CT , £ r and 2ir/(ir — w r ) _1 are com- 
parable to the typical length scale separating KZ de- 
fects (i.e. domain walls). These quantities are given 
by hxG and lnA via Eas. (|H^|) . |STn|l . (|CTT) . Ea. (|73|) 
gives In A = — acosh(l — 2z*/z) which has a known 
large-z expansion. Using the integral representation in 
Eqs.(E3,(E5l,<SSl we obtain 

InG = ln(l - 2e-") - -f Z dx &C0S{1 - X/7TZ) (8.13) 
it Jo e x — 2 

after integrating once by parts and changing variables 
x — > 7rz(l — a;). For large z, we can drop the In term, 
which is exponentially small, expand acos in 1/z, and 
integrate term by term with the upper limit at infinity 
to obtain the large-z expansion. This procedure yields 



- Q (2ttz)" +1 /2 



^ /ln2 

n=0 v 



ri+1/2 



Here the coefficients A n , B n are given by 

r(n+i) 2 RcLi n+3/2 (2) 



A, 



B n 



7T 3 /2r(n+ 1) 

£(rH-l) 

7r 1 /2( n+ |)r(n + l) 



(8.14) 
(8.15) 

(8.16) 
(8.17) 



where Lij,(a;) is the polylogarithm function and T(x) is 
the gamma function. These expressions exhibit the scal- 
ing l a ,Ti (7T — w) _1 cx z 1 / 2 expected from the KZ picture. 
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FIG. 7: Magnetization m z , Ea. l|9.3|l . as a function of the 
inverse sweep speed z/z*. Note the smooth transition from 
small to large z. 



IX. MAGNETIZATION 

The correlators of a x are much simpler to analyze 
since they are composed of a fixed number of Majorana 
fermions: 



u% = A X B X 
&x@x+n = A x B x A x + n B x _i_ r 



(9.1) 
(9.2) 



Averaging these expressions with the help of Ea. (|4.12|l . 
we find 

m z = (a 3 x ) = 1 - 2e- nz I (TTz) (9.3) 

(44+2n) = (4) a -(2e"" In{Kz)f. (9.4) 

where n ^ and, as before, the brackets (...) denote 
the averages in the full lattice. The magnetization m z is 
plotted in Fig. 

The behavior of magnetization at fast and slow sweep 
is given by the small-z and large-z asymptotic: 



-l + 2z + 0(z 2 ) 



z < 1 



1 - (2/ttz) 1 /2 + 0( z -3/2) 2>1 



(9.5) 



where z — ttz. The asymptotic expansion in Eq. l|9.5[l is 
in integer powers at small z, and in negative half- integer 
powers at large z. The small z limit corresponds to fast 
sweep, and so the magnetization deviates little from the 
initial H(t — > — oo) ground state value of m z — — 1. In 
contrast, large z describe slow sweep when the magne- 
tization follows the dynamical field h(t) nearly adiabat- 
ically, and thus m z approaches the H(t — ► +oo) ground 
state value m z = +1. 

The magnetization correlator (cr x a x+2n ) is also a 
smooth function of z. Subtracting (<J X ) 2 , we obtain the 
irreducible (connected) correlator 



r 3_3 \ 

7 x a x+2n/ 



(2 e -^/„(^z)) 2 . (9.6) 
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Correlations of magnetization at distant points are given 
by the large-n expansion of D n at fixed z. We obtain 

e-™/„(7rz) = f e **("*o-i) e ineM = ( 2 ^)-i/2 e -« 2 /2* 

where we used an expansion near the saddle point, 1 — 
cos 8 = \Q 2 + 0(9 4 ), and a gaussian approximation for 
the integral over 9. This gives an asymptotic behavior 

D n = -—e- n *' z (9.7) 

■KZ 

The correlation length, which is very short at small z 
(fast sweep), grows as £ cx z 1 / 2 at large z (slow sweep), 
in agreement with KZ picture. 



X. CONCLUSION 

This article presents an exact solution for a quantum 
spin chain driven through quantum critical points. We 
consider an anisotropic XY chain in a time-dependent 
transverse field hit) that drives the system from a disor- 
dered paramagnetic phase at early times into an ordered 
Ising phase, and back into the paramagnetic phase at 
late times, crossing two quantum critical points along 
the way. We construct an exact many-body evolution 
operator in fermionizcd representation with the help of 
Landau-Zener transition theory, and use it to study the 
evolved state. It is found that the evolved many-body 
state, while technically a pure state, acquires local prop- 
erties of a mixed state. The emerging nonequilibrium 
steady state is characterized by finite entropy density, 
which is a function of the sweep speed. The transforma- 
tion of a pure state into an entropic state, resulting from 
intrinsic decoherence, is analyzed via coarse-graining in 
momentum space. 

Correlation functions in the final entropic state are cal- 
culated using the method of Toeplitz determinants. We 
present exact results for the the asymptotic behavior of 
spin correlators at large spatial separation. The corre- 
lation length dependence on the sweep speed is found 
to be consistent with the Kibble-Zurek —1/2 power law 
scaling. We characterize the crossover behavior in which 
the correlation functions, monotonic at fast speed, ac- 
quire oscillatory spatial dependence at slow speed. The 
critical speed for this transition is found near which the 
the correlation function parameters exhibit nonanalytic 
behavior. 
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APPENDIX A: TOEPLITZ DETERMINANT 
ASYMPTOTICS 

Toeplitz matrices, having constant diagonals, and their 
determinants, Ea. (|6.11[) . arise in many mathematical and 
physical problems. In particular, one is often interested 
in the large-n behavior of D n [f]. Toeplitz determinant 
asymptotics form basis for a number of rigorous results, 
being particularly useful in the computation of various 
quantities in two-dimensional Ising model (see for exam- 
ple Ref. |32J). However, the rather daunting mathemat- 
ical literature on the subject has led to some confusion 
on the status and use of various results such as Szego's 
limit theorem and generalizations of the Fishcr-Hartwig 
conjecture (for example, see Chap. 10 of Ref. [13 )■ We 
give a formulation of Toeplitz determinant asymptotics 
that unifies all previously known results, both conjec- 
tured and mathematically rigorous, and extends them to 
a larger class of Toeplitz determinants. 

The central quantity in the study of Toepltiz determi- 
nants is a function of complex variable, called generating 
function, which is specified for some contour C that en- 
closes the origin once. The generating function fc{Q 
integral over C gives the matrix elements /„ via 

/n=||f"M) (A.l) 

The theory of Toeplitz determinants links the large-n be- 
havior of D n [f] to the analytic structure of fc(Q, in 
particular its singularities. Here we wish to stress two 
points. The first point is the importance of specifying the 
contour C in relating /c(£) to /„ as it gives an explicit 
distinction between singularities inside, outside, and on 
C . This point is well-known in the literature where C is 
taken to be the unit circle and figures prominently in the 
derivation of known results on Toeplitz determinants. 

The second point is the freedom to deform C to C 
when /(£) is analytic between the two contours. This 
point was briefly mentioned in Ref. |33l ] and used to 
obtain the behavior of spin correlation functions of the 
two-dimensional Ising model above the transition tem- 
perature. The freedom to deform C in a general setting 
is a key element of our proposed extension of Toeplitz 
asymptotics. 

We consider the class of generating functions, first 
studied by Fisher and Hartwig [34] , which are given by 

MO = e ff(5) r II (1 - (1 - ^pT 1 )* (A.2) 

p 

with H(£) — h + (t;) + h-(£) + ho, where h + (h ) is an- 
alytic inside and on C (outside and on C) satisfying 
/i_i_(0) = (/i_(oo) = 0), and m is an integer winding 
number. The roots X p are on C and give power-law sin- 
gularities with exponent a p (f3 p ). However, this represen- 
tation of /c(£) is not unique in the sense that different 
choices of C, /iq, h±, \ p , a p , and (3 P will give the same 
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matrix elements /„. For fixed C, the formal identity 

l=(-r 1 A) n (l-A- 1 £) n (l-ATT n (A.3) 
for integer n shows that the transformation 

p 

m — » m — n p (A. 5) 

p 

a p — > a p + n p (A. 6) 

/3 p -> /3 p -n p (A.7) 

gives a generating function with the same Fourier coef- 
ficients / n , Ea. l|A.l|) . as those obtained for the original 
function fc(£,)- Under such a transformation, while the 
parameters ho, a p , (3 P , m change, the Toeplitz matrix is 
preserved. The consequences of such transformation was 
first pointed out by Basor and Tracy in |35|. who noted 
that all different generating function representations con- 
tribute to the asymptotics. 

Let us now consider deformations of the contour C. We 
note that Eq. (|A.2fl allows singularities to be on C. Since 



the matrix elements f n given by Eq. ljA.ljl must remain 
the same upon deforming C to C", such a deformation 
must not enclose any singularities, but C can possibly 
pass through additional singularities that C does not. In 
the representation ljA.21) , singularities strictly outside (in- 
side) of C are described by h_ and h+ while singularities 
on C are described by the roots A p . By appropriately 
deforming C to C , we can move power-law singularities 
from h + and h- and include them in additional roots X' p 
on C 

The most general result in the literature is for C fixed 
to be the unit circle but taking into account the trans- 
formations of Eqs. lA.4jl - l)A.7|) . It was first proposed by 
Basor and Tracy |35| and is known as the generalized 
Fishcr-Hartwig conjecture. Each representation given by 
Eqs. l|A.4f) - i|A.7|) gives a contribution to D n [f] of the form 



5 mt0 Ae h ° n rfip a pPp (A.8) 

with the prefactor 



A = cxp ( khkh- 



n 



G(l + a p )G(l + P P ) 
G(l + a p + (3 P ) 



P h-{X V )-I3 p h + {X V ) -Q L_ V N 



(A.9) 



p 



where G(x) is the Barnes G-function j^a] which satisfies 
G(x + 1) = T{x)G{x) and 

k=l 

The constraint 5 m> o in Ea. (|A.8ll means that the contri- 
butions for non-zero winding numbers m are not of the 
above form but decay faster than nJ 1 for all real rj < 0. 
The asymptotic of D n [f] is obtained by summing the 
terms which give the leading contribution for large n. 

This conjecture has been proven rigorously in some 
cases. The case for arbitrary a p and (3 P , but such that 
only one representation contributes to the lea ding term, 
has only been proven relatively recently j3|| l37j . The 
case for positive integer a p and f3 p but with multiple rep- 
resentations contributing at leading order was proven by 
Bottcher and Silbermann |3l| . 

The generalized Fisher-Hartwig conjecture as stated 
above gives the asymptotics of D n [/] as a sum of contri- 
butions from each equivalent representation of the gen- 
erating function fc(Q; but with C fixed to be the unit 
circle. The natural extension is to also allow arbitrary de- 
formations of C to C that may touch but not cross the 
singularities, and then sum over the leading contributions 
from the additional equivalent representations. This pro- 



cedure can be concisely expressed as follows. One writes 
down the generating function in the form of Ea. l|A.2|) 
with the power-law singularities for X p arbitrarily dis- 
tributed in the complex plane. Then one generates the 
equivalent representations via Eas. ijA.4(l - l|A.7|l and sums 
over the contributions given by Ea. ljA.8|l . In practice 
only the singularities closest to the unit circle need to 
be considered. This is because under the transforma- 
tion of Eas. ljA.4(l - l|A.7|) . e h ° gets multiplied by powers of 
A p . Since the contribution to D n [f] given by Ea. i|A.8|) 
contains e h ° n , the singularities far from the unit circle 
generically give subleading contributions. 

This proposed extension of the generalized Fisher- 
Hartwig conjecture is particularly useful for generating 
functions with power-law singularities off of the unit 
circle but non-zero winding number. The generalized 
Fishcr-Hartwig conjecture is not applicable in this case 
due to the presence of winding numbers. [Note 5 m Q in 
Ea. HA.8jl .] However, by using the freedom to deform 
the contour C to the singularities off of the unit cir- 
cle we can absorb the winding number into a p , (3 P via 
Eas. l|A.4f) - i|A.7|) . After that the zero winding number re- 
sult, Ea. HA.8p . can be used. 

Essentially, the generalized Fisher-Hartwig conjecture 
relates power-law singularities on the unit circle to the 
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asymptotic behavior of Toeplitz determinants. Our pro- 
posed extension just states that the same relation holds 
for power-law singularities with a generic location in the 
complex plane. The literature on Toeplitz determinants 
mostly considers singularities on the unit circle, with an 
exception of the result obtained by Day [2S| for rational 
generating functions 



(A.10) 



function is clearly of the form of Ea. ljA.2j) . The corre- 
sponding Toeplitz determinant can be evaluated explic- 
itly and the result is given exactly by generating all equiv- 
alent representations using all the roots via Eas. (jA.4j) - 
(|A.7() and summing over the corresponding contributions 
of Ea. HA.8j l. This provides evidence that the extension 
proposed here holds in a general setting, although we ex- 
pect it to give only the leading asymptotic contribution 
and not the exact determinant in this case. 



where X p and p q are arbitrary in the complex plane. This 
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